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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC)  UNITS  OP  MEASUREMENT 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI  (metric)  units  as  follows: 


Multiply 

By 

To  Obtain 

feet 

0.3048 

metres 

inches 

2.54 

centimetres 

inches 

0.0254 

metres 

pounds 

4.4822 

newtons 

tons 

8.896 

kilonewtons 

pounds  j>er  square  foot 

47.8803 

pascals 

pounds  per  square  foot 

0.04788 

kilopascals 

Founds  i>er  square  inch 

6.8948 

kilopascals 

tons  per  sqaure  foot 

95.76 

kilopascals 

tons  per  sqaure  foot 

0.976 

kg/cm^ 
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SOLUTION  OF  SOIL-STRUCTURE  INTERACTION  PROBLEMS 
BY  COUPLED  BOUNDARY  ELEMENT-FINITE  ELEMENT  METHOD 


Introduction 

Soil-Structure  interaction  problems  in  geotechnical  engineering  involve  the  solution  of  boundary 
value  problems  consisting  of  two  domains,  a  finite  domain  representing  the  structure  and  a  semi- 
infmite  domain  representing  the  soil  into  which  the  structure  is  embedded.  These  two  distinct  domains 
are  separated  by  a  boundary  called  the  soil-structure  interface.  Figure  1  illustrates  a  typical  soil 
structure  interaction  problem.  The  aim  of  solving  the  boundary  value  problem  is  to  determine  the 
displacements  and  stresses  in  both  domains.  In  general,  it  is  very  difficult  to  obtain  a  closed  form 
solution  to  the  above  problem  and  in  most  cases  the  solution  must  be  obtained  numerically.  Recently, 
the  finite  element  method  has  been  widely  applied  to  the  solution  of  this  class  of  boundary  value 
problems.  The  main  advantage  of  this  method  is  its  capability  to  solve  nonlinear  boundary  value 
problems  with  arbitrary  geometry,  boundary  conditions  and  constitutive  properties.  One  limitation  of 
this  method  however  is  that  it  can  only  solve  finite  domain  problems.  Thus,  infinite  domain  problems 
are  approximately  solved  by  taking  a  considerably  large  domain  and  assuming  the  displacements 
vanish  along  the  far-field  boundary  distant  from  the  region  of  interest.  More  recently,  the  boundary 
element  method  has  been  used  to  solve  such  infinite  domain  problems.  For  linear  problems,  the 
boundary  element  method  results  in  a  smaller  system  of  equations  with  much  greater  accuracy 
particularly  in  infinite  and  semi-infinite  domain  problems.  One  limitation  of  this  method  is  that 
problems  involving  non-hcmogencous  and  non-linear  domains  result  in  a  more  complicated  solution 
procedure.  A  compromise  between  the  two  methods  would  be  to  model  the  non-linear  and  non- 
homogeneous  domain  in  the  near  field  region  of  interest  using  finite  elements  and  to  model  the  far-field 
boundary  using  boundary  elements  assuming  that  the  distant  domain  is  linear  and  homogeneous.  The 
resulting  problem  is  a  coupled  model  involving  both  boundary  element  and  finite  element  systems.  The 
main  setback  in  this  approach  is  that  the  symmetry  and  bandedness  of  the  system  of  equations 
resulting  from  the  finite  element  method  is  destroyed  by  the  boundary  element  method.  Thus,  the 
matrix  resulting  from  the  coupling  of  both  systems  is  unsymmetric  as  well  as  fully  populated.  The 
purpose  of  this  paper  is  to  show  that  use  of  boundary  elements  for  modeling  the  far  field  in  the  finite 
element  analysis  results  in  a  significant  improvement  in  accuracy  and  such  a  technique  can  be 
efficiently  implemented. 
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Figure  1  A  Typical  Soil-Structure  Interaction  Problem 


Boundary  Element  Method 


The  direct  boundary  element  method  is  a  numerical  technique  for  solving  boundary  value 
problems.  In  this  method,  the  governing  differential  equations  are  converted  into  an  equivalent  integral 
statement  by  successively  applying  the  divergence  theorem.  This  integral  statement,  called  the 
boundary  integral  equation,  for  an  elastostatic  boundary  value  problem  is  given  by  (Brebbia  [1984]) 

Cij(s)ii_,(s)+  /  p’j{s,q)  Uj(q)dr  =  f  n*j(s,q)  Pjiq)dr  (1) 

Jr  dp 

The  displacements  and  tractions  Uj(q)  and  Pj{q),  are  the  jth  component  at  p>oint  q  on  the  boundary 
r.  The  displacements  and  tractions  and  p’j(s,9),  are  the  fundamental  solution  corresponding 

to  the  jth  component  at  point  q  due  to  a  unit  force  applied  at  ptoint  s  in  the  >th  direction  where  both 
{x>int8  s  and  q  are  on  the  boundary  P  and  C',j(s)  is  a  2  x  2  matrix  which  is  a  function  of  the  pwint  s. 

The  basic  steps  in  the  boundary  element  method  involve: 

(a)  The  boundary  P  is  discretized  into  a  series  of  elements  of  which  displacement  and  tractions  are 
chosen  to  be  piecewise  interpolated  between  the  element  nodal  points; 

(b)  The  boundary  integral  equation  is  applied  in  discretized  form  to  each  nodal  point  s  on  the 
boundary  P  and  the  integrals  are  computed  (usually  by  a  numerical  quadrature  scheme)  over 
each  boundary  element.  A  system  of  N  linear  algebraic  equations  involving  the  set  of  N  nodal 
tractions  and  N  nodal  displacements  is  therefore  obtained; 

(c)  Boundary  conditions  are  imposed  and  consequently  N  nodal  values  (traction  or  displacement  in 
each  direction  p>er  node)  are  prescribed.  The  system  of  N  equations  can  therefore  be  solved  by 
standard  methods  to  obtain  the  remaining  boundary  conditions. 

In  what  follows,  each  of  the  above  stepra  shall  be  explained  more  thoroughly  together  with  an 
application  to  the  solution  of  two  dimensional  elastostatic  problems  using  linear  elements. 

For  the  discretization  of  the  boundary  integral  equation,  the  boundary  P  is  approximated  by 
>ing  a  series  of  elements.  The  Cartesian  coordinates  {x,y)  of  pwints  located  within  each  element  P- 
ire  expressed  in  terms  of  the  lagrangian  interpolating  functions  rpj  and  the  nodal  coordinates  (xj,  yj)  of 
the  element  by  the  following  relation 
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(2) 


j  =  1  ]  =  \ 

Similarly,  the  boundary  displacements  and  tractions  are  approximated  over  each  element  in  terms 
of  the  lagrangian  interpolating  function  rp-  in  the  following  relation 

=  Z!  ^ij^j  Pi=Yl  (3) 

i  =  1  j  =  1 

where  u,-  and  p-  is  the  ith  component  of  the  displacements  and  tractions  along  the  element  and  t/,j  and 
P^j  is  the  ith  component  of  the  nodal  tractions. 

Assuming  that  the  boundary  P  is  discretized  into  L  elements  P^,  the  substitution  of  the  above 
approximations  into  eqn.  (1)  results  in  matrix  form 

C(sM^i)+  J  ^p’(s-,q)^dr^Uiq) 


u’(Si,q)xl>dryiq) 


(4) 


for  a  boundary  node  s-  . 


In  general,  the  integrals  presented  in  the  above  equation  must  be  computed  numerically. 
However,  for  the  case  when  the  singular  node  i  lies  within  the  element  P j,  the  integral  becomes 
improper  because  of  the  singularity  in  both  p*  and  w*.  The  special  treatment  of  theses  integrals  shall 
be  described  later.  For  the  case  when  the  singular  pwints  is  outside  the  element  P^  the  integral  can  be 
readily  replaced  by 


/ 

I 


^dP  = 


u  fPdP  = 


n  ^ 

J  k  =  l 


r\  ^ 

I  u"  Ji>\j\dr]  = 

-1  k  =  1 


(5) 

(6) 


where  K  is  the  number  of  Gaussian  integration  points,  is  the  weighting  factor  associated  with 
them,  T)  is  the  local  coordinate  system  and  I J  I  is  the  Jacobian  of  the  transformation  given  by 
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(7) 


Application  of  eqn.  (4)  to  all  nodes  results  in  a  system  of  equations  of  the  form 


(C  +  H)V  =  GP 
or 

HU  =  GP 


(8) 

(9) 


where  H  —  C  +  H  and  C  is  a  diagonal  matrix  consisting  of  the  constants  C(s,)  mentioned  previously. 

For  the  sptecial  case  when  the  p>oint  s,-  is  one  of  the  nodal  points  of  the  element  P^  ,  the  right- 
handside  term  involving  u*  must  either  be  integrated  analytically  as  is  p)ossible  for  constant  and  linear 
elements  (Brebbia  [1978]  and  [1984],  Crouch  [1983]),  or  for  higher  order  elements,  the  logarithmic  part 
must  be  separately  integrated  using  a  special  numerical  quadrature  (Cristescu  [1978]).  Other  methods 
involve  developing  techniques  where  the  singular  terms  cancel  (Zeng  [1990]). 

It  should  be  noted  that  for  the  total  computation  of  H,  the  leading  diagonal  submatrices  which 
correspond  to  the  term  C,j,  the  principal  value  of  the  integral  in  eqn.  (4)  can  be  calculated  by 
imposing  to  the  matrix  equation  the  condition  that  rigid  body  translations  esult  in  zero  tractions. 
Thus  result  of  this  is  that  for  eqn.  (9),  H  can  be  computed  using  the  following  relations  (Brebbia 
[1978]) 

NN 

^aa  =  -  E  (10) 

q  ^  a 

for  finite  bodies  and 

NN 

(H) 

q  ^  a 

for  infinite  and  semi-infinite  domains.  In  the  above  expressions  H-  refers  to  a  2x2  submatrix 
corresponding  to  the  evaluation  of  the  fundamental  solution  at  nodes  i  and  j  and  I  is  a  2  x  2  identity 
matrix. 

Fundamental  Solution 

For  a  general  two  dimensional  plane  strain  elastostatic  boundary  value  problems,  the 
fundamental  solution  cissociated  with  a  unit  load  applied  at  a  point  s  within  infinite  plane,  the 
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displacements  and  tractions  at  point  q  are  given  by  (Timoshenko  [1970]) 


8ir(]  —  i/)G 
4ir(l 


where  the  subscripts  i  and  j  refer  to  the  directions  of  the  applied  unit  force  and  the  resulting 
component  of  the  corresponding  quantity  respectively.  Also, 
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•j 


0  for  i  ^  j 
1  for  t  =  j 


and  where  r  is  the  distance  between  points  s  and  g,  r,-  is  the  ith  component  of  r  and  n  is  the  unit 
vector  normal  to  the  boundary  jT  at  point  q,  and  1/  is  the  Poisson’s  ratio.  In  most  boundary  element 
literature,  this  solution  if  referred  to  as  the  Kelvin  solution.  The  above  fundamental  solution  is  suitable 
for  solving  both  finite  and  infinite  domain  problems.  For  soil-structure  interaction  problems  however, 
there  is  a  need  to  model  the  far  field  as  an  infinite  half-plane.  The  fundamental  solution  corresponding 
to  this  problem  was  obtained  by  Melan  [1932].  However,  in  the  original  work,  there  seem  to  be  errors 
which  were  pointed  out  by  Telles  [1981].  The  analytical  solution  to  this  problem  consists  of  superposing 
a  complementary  solution  (  )*  on  the  Kelvin  solution  (  .  Thus  the  solution  can  be  expressed  as 


u‘  =  u;  +  u; 


P*  =  pI  +  p: 


The  complete  and  correct  expression  for  the  complementary  solution  is  given  in  Appendix  B.  In  the 
proposed  solution  procedure,  the  Melan  fundamental  solution  is  used  to  account  for  the  semi-infinite  far 
field  domain.  This  aspect  of  the  procedure  contrasts  with  other  researchers,  Vallabhan  [1986]  and 
Georgiou  [1981]  for  example,  who  use  the  Kelvin  fundamental  solution  in  their  boundary  element 
procedure  to  model  the  far  field.  It  should  be  pointed  out  that  the  result  of  adopting  the  Kelvin 
solution  is  that  the  far  field  domain  must  be  treated  as  a  large  finite  body  and  the  displacements  must 
be  assumed  to  vanish  as  points  away  from  the  region  of  interest  are  considered. 


Linear  Boundary  Elements 


Among  the  different  types  of  elements  that  can  be  employed  in  the  numerical  disretization  of  the 
boundary  integral  equation,  the  linear  element  has  been  found  to  give  acceptable  accuracy  without 
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r 


requiring  much  computer  effort  for  the  solution  of  the  examples  presented  here.  The  geometry  of  the 
element  is  represented  by  a  straight  line  of  length  /  and  |  J|  =  //2.  Also  the  interpolating  functions  are 
given  by 


= 


(1-0 

2 


(12) 


It  can  be  seen  that  the  discretization  of  the  boundary  integral  equation  gives  rise  to  two  element 
matrices  h  and  g  (corresponding  to  B  and  G  respiectively)  of  order  2  x  4  in  the  following  form 


I 


(13) 

(14) 


For  the  sp>ecial  case  where  the  singular  node  is  coincident  with  one  of  the  end  nodes  of  the  element,  the 
coefficients  of  the  g  matrix  can  be  computed  analytically  to  avoid  significant  error  introduced  by 
numerical  integration  of  the  improper  integral.  For  the  Kelvin  solution,  the  following  expression  is 
obtained 

=  (i6x(i-i/)c;] 


for  the  half-plane  solution,  the  complementary  part  is  given  by 


{2(1  -..)/  -h  0.5  -ln(/)Ii, .  -h  lili +  a(l- 6, -y  9  (1  - 


where  the  subscripts  i  and  j  indicate  the  p>osition  of  the  coefficient  in  the  2x2  submatrix  k  of  g  and  n 
indicates  the  singular  node.  Also 

_  f  0  for  »  j 
*^“jlfor.  =  j 

{-1  for  j  =  1 
1  for  j  =  2 

9  =  arctan(/2//i)  (  -  ^/2  <9  <  x/2) 
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and  where  I-  is  the  ithe  component  of  the  element  length. 


Coupling  the  Boundary-Finite  Element  System 

For  an  elastostatic  boundary  value  problem,  the  system  of  linear  equations  resulting  from  the 
finite  element  method  can  be  written  in  the  form 

KU  =  F  (15) 

where  K  is  the  stiffness  matrix,  u  is  the  vector  of  nodal  displacements,  and  F  is  the  vector  of  nodal 
forces.  The  resulting  system  of  equations  resulting  from  the  boundary  element  method  can  be  expressed 
in  the  form 

SU  =  GP  (16) 

where  H  and  G  are  the  matrices  of  influence  coefficients,  «  is  the  vector  of  nodal  displacements,  Md  p 
is  the  vector  of  nodal  tractions  in  the  boundary  element  system.  The  nodal  traction  can  be  converted 
into  equivalent  nodal  forces  using  an  appropriate  definition  of  the  form 

=  (17) 

where  IV  is  a  square  transformation  matrix  which  converts  the  nodal  tractions  to  equivalent  nodal 
loads.  Using  the  above  definitions,  the  modified  boundary  element  system  is  of  the  form 

HU  =  GN-^F  (18) 

which  can  be  rewritten  as 

KU  =F  (19) 

where 

K  =  G-^NH  (20) 

To  describe  the  procedure,  assume  that  the  finite  element  system  can  be  partitioned  into  those 
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nodes  which  are  included  only  in  the  finite  element  system,  and  those  nodes  common  to  both  the  finite 
element  and  boundary  element  system.  The  partitioned  system  of  equations  thus  can  be  written  as 


A'// 

^6/  ^66 


Ut 


(21) 


where  and  /’j,  are  the  displacements  and  nodal  forces  respectively  lying  on  the  interface  of  the  finite 
element  and  boundary  element  domains  while  tiy  and  are  the  displacements  and  nodal  forces 
respiectively  not  lying  on  the  finite  element  and  boundary  element  interface. 

The  boundary  element  system  can  be  written  as 


KU  =  F 


(22) 


The  variables  in  the  two  systems  are  related  by  continuity  and  equilibrium  which  require  that 


V^=U  (23) 

Ffc  +  F  =  0  (24) 


or 


F=-F, 


(25) 


Thus,  applying  the  conditions  mentioned  in  Eqns. 


(9)  and  (10)  yield  the  coupled  system  of  the  form 


K 

K 


Jf 

6/ 


^/6 

Kbb  +  F 


(26) 


It  should  be  noted  that  the  stiffness  matrix  K  is  symmetric  positive  definite  and  banded  while  the 
stin'ness  matrix  for  the  modified  boundary  element  method  K  is  non-symmetric  and  fully  populated. 
Thus,  because  of  the  boundary  element  coupling,  the  symmetry  and  bandedness  of  the  stiffness  matrix 
in  the  finite  element  method  is  lost. 

Some  researchers  (Vallabhan  [1986],  Georgiou  [1981])  have  tried  to  retain  the  symmetry  of  the 
stiffness  matrix  by  developing  an  approximate  stiffness  matrix  for  the  modified  boundary  element 
method.  One  way  of  doing  this  is  by  taking  only  the  symmetric  part  of  the  matrix  and  discard  the 
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skew-symmetric  part.  Thus  the  approximate  stiffness  matrix  is  given  by 


(27) 


where  k-  is  the  (»,  j)th  element  of  the  K  matrix. 

In  the  process  of  this  research,  the  effect  of  discarding  the  skew-symmetric  portion  of  the  H  and  G 
matrix  on  the  accuracy  of  the  solution  was  investigated.  This  consisted  of  solving  typical  elasticity 
problems  using  both  the  unsymmetric  and  symmetric  H  and  G  matrices  in  the  boundary  element 
solution.  Results  from  this  investigation  indicate  that  if  the  Kelvin  solution  is  used  and  the  problem  is 
solved  without  taking  advantage  of  the  geometric  symmetry  of  the  problem,  there  is  no  significant 
effect  due  to  discarding  the  skew-symmetric  part.  However,  if  geometric  symmetry  is  taken  into 
account  by  reflecting  the  singular  nodes,  significant  errors  occur.  If  the  Melan  solution  is  used,  it  was 
determined  that  the  error  due  to  discarding  the  skew-symmetric  psirt  is  dependent  upon  the  value 
assigned  to  poisson’s  ratio  u.  That  is,  as  i/  approaches  0.5,  the  error  increases  regardless  of  whether  or 
not  symmetry  is  taken  into  account.  It  should  be  noted  however  that  when  the  boundary  element 
method  is  coupled  with  the  finite  element  method  in  the  procedure  described  above,  the  effects  of 
discarding  the  skew-symmetric  part  on  the  accuracy  of  the  final  results  are  minimized.  There  are  two 
reasons  for  this.  First,  because  the  elements  of  the  boundary  element  stiffness  matrix  derived  in  eqn. 
(20)  are  considerably  smaller  that  the  finite  element  stiffness  matrix  into  which  they  are  assembled  and 
thus  the  symmetry  of  the  finite  element  matrix  decreases  the  ratio  of  the  skew-symmetric  psirt  over  the 
symmetric  part  of  the  overall  stiffness  matrix.  Secondly,  it  was  observed  that  because  the 
transformation  matrix  N  is  symmetric,  the  resulting  stiffness  matrix  K  is  more  symmetric  than  that 
corresponding  to  the  boundary  element  method  alone  and  given  by 

K  =G~^H  (28) 

The  coupling  procedure  can  be  summarized  in  to  the  following  steps: 

STEP  I:  From  the  problems  geometry  and  material  parameters  given  form  the  H,  G  and  N  matrices. 
The  formulation  of  the  N  matrix  is  given  in  Appendix  C. 

STEP  2:  Perform  the  matrix  operations  specified  in  eqn.  (20)  to  obtain  the  boundary  element  stiffness 
matrix  K. 

STEP  3:  If  only  the  symmetric  part  of  K  is  to  be  used,  discard  the  skew-symmetric  p)art  using  the 
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procedure  described  in  eqn.  (27). 


STEP  4:  Assemble  the  stifTness  matrix  A'  into  the  global  stiffness  matrix. 

STEP  5:  Solve  the  resulting  system  of  equations  as  would  be  done  for  a  purely  finite  element  system. 
However,  if  the  full  K  is  used,  an  usymmetric  solution  technique  must  be  used. 

STEP  6:  Using  the  displacements,  compute  the  stresses  as  would  be  done  in  a  purely  finite  element 
system. 

It  should  be  noted  from  the  above  discussion  that  once  the  equivalent  stiffness  matrix  for  the 
boundary  element  system  is  obtained  it  can  be  assembled  into  the  finite  element  system  in  the  same 
way  finite  elements  are  assembled.  Thus,  the  coupling  procedure  proposed  can  be  viewed  as  treating  the 
boundary  element  system  as  a  substructure  and  obtaining  the  stifTness  matrix  corresponding  to  this 
substructure.  The  substructure  is  treated  as  a  special  finite  element  which  can  be  incorporated  into  the 
element  library  of  an  exiting  finite  element  program.  In  order  to  obtain  a  system  of  equations  where 
the  degrees  of  freedom  common  to  both  the  finite  element  and  boundary  element  system  occur  as 
described  in  eqn.  (26),  a  procedure  of  renumbering  the  equations  is  used.  The  process  of  renumbering  is 
described  in  greater  detail  in  the  following  section.  The  main  purpose  in  separating  these  degrees  of 
freedom  is  to  minimize  the  band  width  that  comes  from  the  dense  nature  of  the  substructure’s  stiffness 
matrix. 

Program  Implementation 

A  coupled  boundary  element-finite  element  program  BEFEC  was  developed  to  verify  the  validity 
of  the  proposed  coupling  technique.  The  finite  element  portion  of  BEFEC  was  derived  from  static 
portion  of  DLEARN,  a  program  developed  by  Hughes  (1988).  The  boundary  element  portion  was 
derived  from  the  elastostatic  2D  program  developed  by  Brebbia  [1984]  and  contains  additional  features 
such  as  the  half-plane  fundamental  solution.  The  system  of  equations  are  assembled  into  an  active 
column  form  where  only  the  non-zero  quantities  in  the  global  stiffness  matrix  occurring  beneath  the 
skyline  are  stored.  The  scheme  of  storage  is  one  of  the  most  efficient  in  terms  of  amount  of  memory 
used  as  well  as  computer  time  required.  The  only  modification  in  the  algorithm  used  to  assign  the 
equations  numbers  is  that  two  passes  instead  of  one  are  used.  In  the  first  pass,  equation  numbers  are 
assigned  to  those  degrees  of  freedom  corresponding  to  nodes  which  are  not  on  the  far  field  boundary.  In 
the  second  pass,  the  nodes  corresponding  to  the  far  field  boundary  are  assigned  equation  numbers.  This 
algorithm  assures  that  the  fully  populated  boundary  element  stiffness  matrix  is  confined  to  the  last 
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equation  and  prevents  these  terms  from  dispersing  throughout  the  entire  stiffness  matrix.  The  result  is 
that  the  column  heights  are  minimized.  The  program  also  has  the  option  to  discard  the  skew- 
symmetric  po.'iion  of  the  boundary  element  stiffness  matrix  thus  resulting  in  an  approximate 
symmetric  boundary  element  stiffness  matrix  as  done  in  eqn.  (27).  Vallabhan  [1986]  suggests  the 
following  condition  be  satisfied  before  discarding  the  skew  symmetric  part  in  order  that  the  resulting 
error  not  be  significant 

maximum  |  |  ^  (*'j) 

where  o-  are  the  elements  of  the  global  stiffness  matrix.  It  was  noted  that  in  the  problems  considered, 
the  above  norm  was  much  less  than  1%.  Vallabhan  [1986]  states  that  the  above  norm  becomes  smaller 
as  the  size  of  the  problem  increases. 

It  should  be  noted  that  in  the  solution  of  soil-structure  interaction  problems,  while  similar  in 
many  ways,  the  above  solution  procedure  has  three  major  differences  from  previously  developed 
coupling  techniques  to  solve  such  problems  (Vallabhan  [1986],  Georgiou  [196iJ).  First,  the  far  field 
domain  is  modeled  using  boundary  elements  based  on  the  half-plane  (Melan)  solution.  This  approach 
has  be  commonly  used  in  dynamic  soil  structure  interactions  problems  (Manolis[1988],  Wolf[1988]). 
However,  this  is  not  so  in  the  static  case  where  the  far  field  domain  is  modeled  as  a  large  finite  body 
using  boundary  elements  based  on  the  infinite  plane  (Kelvin)  solution.  Also,  in  the  works  cited, 
constant  boundary  elements  are  used  as  compared  to  the  linear  boundary  element  used  in  the  current 
formulation.  Lastly,  in  the  works  cited,  no  renumbering  is  performed  on  the  system  of  equations, 
making  it  necessary  to  store  the  full  system  of  equations  making  the  solution  procedure  costly  in  terms 
of  memory  required.  It  should  also  be  stressed  that  many  existing  finite  element  programs  use  the 
active  column  storage  scheme.  Thus  the  proposed  coupling  technique  can  easily  be  integrated  into  such 
program  since  it  merely  involves  adding  a  new  element  to  the  existing  element  library  of  these 
programs. 

Program  Validation 

To  validate  the  proposed  solution  method,  two  linear  elastic  half-plane  problems  were  solved.  For 
each  problem,  five  different  analyses  were  performed.  This  involved  first  conducting  a  finite  element 
analysis  using  a  suitable  mesh  and  assuming  the  far  field  displacements  vanish.  A  second  finite  element 
analysis  was  conducted  using  an  expanded  domain  while  retaining  the  same  refinement  as  the  first 
analysis.  A  third  finite  element  analysis  was  conducted  using  the  same  domain  as  the  first  analysis  but 
more  refined  mesh.  Two  coupled  boundary  element-finite  element  analyses  were  conducted  using  the 
same  mesh  as  in  the  first  finite  element  analysis  but  modeling  the  far  field  effects  using  boundary 
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elements.  In  the  first  coupled  analysis  the  unsymmetric  part  of  the  stiffness  matrix  resulting  from  the 
boundary  element  method  was  considered  while  the  analysis  neglected  this  part  of  the  stiffness  matrix. 


Problem  1:  This  problem  consists  of  a  concentrated  load  applied  to  a  half  space.  To  solve  this 
problem,  four  different  analyses  were  considered  namely:  (1)  a  finite  element  analysis  consisting  of  a  30 
m.  X  3.0  m.  region  using  100  nodes  and  81  elements,  FEM  (100)  [fig.  2]  ,  (2)  a  finite  element  analysis 
consisting  of  a  60  m.  x  60  m.  region  using  169  nodes  and  144  finite  elements,  FEM  (169A)  [fig.  3],  (3) 
a  finite  element  analysis  consisting  of  a  30  m.  x  30  m.  region  using  169  nodes  and  144  finite  elements, 
FEM  (169B)  [fig.  4],  and  (4)  a  coupled  boundary  element/finite  element  analysis  consisting  of  a  30 
m.  X  30  m.  region  with  100  nodes  and  81  elements  together  with  18  boundary  elements  (BEFEC- 
18/100)  and  a  symmetric  solution  (BEFEC2-18/100)  where  the  skew-symmetric  portion  of  the  stiffness 
matrix  was  discarded.  Symmetry  of  geometry  was  considered.  The  material  properties  were  £=3.0 
MPa,  1/  =  0.25.  The  load  applied  had  a  magnitude  of  0.5  KN.  The  complete  analytical  solution  is  given 
is  Appendix  A.  In  comparing  the  results  of  the  numerical  analysis  with  the  analytical  solution,  the 
following  observations  are  made. 

1.  Vertical  Displacements  Along  the  Centerline  Uj^(i  =  0,y):  The  initial  finite  element  analysis 
FEM  (100)  shows  a  significant  underprediction  of  the  displacements  (see  fig.  5).  No  significant 
improvement  is  obtained  by  refining  the  mesh  as  can  be  seen  from  the  results  of  FEM  (169B).  A 
substantial  improvement  is  obtained  when  the  domain  of  the  problem  is  increased  as  shown  by  the 
results  of  FEM  (169A).  Results  of  the  proposed  technique  show  a  very  good  agreement  with  the 
analytical  solution.  It  is  also  observed  that  there  is  practically  no  difference  between  the  unsymmetric 
and  symmetric  solution. 


2.  Horizontal  Displacements  Along  the  Surface  u^{x,y  =  0):  Evaluating  the  analytical  solution,  it 
can  be  seen  that  the  horizontal  displacements  along  the  surface  are  equal  to  a  constant. 


u 
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-Py{l-2U) 

AG 


Because  the  boundary  conditions  are  assumed  to  vanish  away  from  the  load  area,  all  three  finite 
element  solutions  give  answers  which  diverge  from  the  analytical  solution  while  both  the  coupled 
solutions  are  give  much  better  answers  which  more  closely  agree  with  the  analytical  solution  as  can  be 
seen  in  fig,  (6)  . 

3.  Vertical  displacements  along  the  surface.  Results  presented  in  fig.  (7)  show  that  like  the 
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Fig.  2  FEM  (100)  Mesh 
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Fig.  4  FEM  (169B)  Mesh 
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Fig.  5  Problem  1 :  Vertical  Displacements  Along  Centerline 
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X  Coordinate  (m.) 

Fig.  6  Problem  1 :  Horizontal  Displacements  Along  Surface 


displacements  along  the  centerline,  the  conventional  method  of  finite  element  analysis  FEM  (100) 
results  is  a  substantial  disparity  between  the  analytical  and  numerical  solution.  This  disparity  is  little 
affected  by  refining  the  mesh  as  was  done  in  FEM  (169B).  The  prediction  of  vertical  displacements  can 
however  be  improved  by  coupling  the  solution  with  the  boundary  element  method  BEFEC  (100/18) 
and  BEFEC2  (100/18),  or  increasing  the  domain  of  the  problem  as  in  FEM  (169A).  In  this  particular 
case,  the  displacements  obtained  using  the  former  procedure  are  much  more  accurate  than  those  using 
the  coupled  solutir:^.  Almost  no  difference  was  noticed  between  the  symmetric  and  unsymmetric 
solutions. 


4.  Horizontal  Normal  Stresses  (<t^^)  along  the  diagonal  elements.  From  fig  (8)  it  can  be  seen  that 
the  stresses  predicted  by  the  finite  element  method  alone  FEM  (100)  and  FEM  (169 A)  are  substantially 
underpredicted.  It  also  can  be  seen  that  increasing  the  mesh  size  decreases  the  amount  by  which  the 
stresses  are  underpredicted  as  is  seen  by  the  difference  between  two  analyses.  In  comparison,  except  for 
points  near  the  loaded  area  where  the  effects  of  stress  concentration  are  dominant,  both  coupled 
solutions  agree  very  well  with  the  analytical  solution.  To  further  study  the  error  involved  in  the  stresses 
computed  using  the  different  methods  of  analysis,  the  relative  error  was  ploted  versus  the  x-coordinate 
in  fig.  (9).  the  relative  error  is  defined  as  follows 


relative  error  (%)  = 


X  100% 


(30) 


where  «r“  is  the  analytical  solution  and  <r"  is  the  numerical  solution.  It  can  be  seen  that  except  for  the 
loaded  area  (x  <  2.5  m.),  the  coupling  method  yields  stresses  which  are  well  within  2%  of  the  analytical 
solution  while  the  finite  element  analysis  alone  yields  errors  of  as  much  as  35%  in  both  zmalyses 
indicating  that  the  prediction  is  not  improved  by  expanding  the  mesh  size. 


5.  Shear  Stresses  *long  the  diagonal  elements.  Results  of  the  analyses  pierformed  are 

presented  in  fig.  (10).  From  the  graph,  it  can  be  seen  that  the  stress  prediction  by  finite  element 
method  alone  is  reasonably  accurate.  Expanding  the  mesh  as  in  FEM  (169A)  or  coupling  with  the 
boundary  element  method  as  in  BEFEC  (100/18)  and  BEFEC2  (100/18)  have  no  visible  effect  on  the 
resulting  solution  indicating  that  these  methods  are  convergent  with  respect  to  the  accuracy  of  the 
solution.  However,  upon  comparing  with  the  relative  error  as  presented  in  fig.  (11),  it  can  be  seen  that 
the  far  field  solution  is  grossly  in  error.  In  particular,  an  underprediction  of  68%  occurs  in  the  initial 
finite  element  solution  FEM  (100).  This  error  is  substantially  reduced  to  20%  by  expanding  the  mesh 
to  twice  the  initial  size  as  in  FEM  (169A).  The  results  obtained  using  both  the  unsymmetric  BEFEC 
(100/18)  and  symmetric  BEFEC2  (100/18)  coupled  solution  yield  identical  values  within  2%  of  the 
analytical  solution. 
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Fig.  7  Problem  1 :  Vertical  Displacements  Along  Surface 
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Fig.  8  Problem  1 :  Sigma-XX  Along  Diagonal  Elements 
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X  Coordinate  (m.) 

Fig.  9  Problem  1 .  Error  in  Sigma-XX  Along  Diagonal  Elements 
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X  Coordinate  (m.) 

Fig.  10  Problem  1:  Sigma-XY  Along  Diagonal  Elements 
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Fig.  1 1  Problem  1 :  Error  in  Sigma-XY  Along  Diagonal  Elements 


6.  Vertical  Normal  Stresses  {<^yy)  along  the  diagonal  elements.  Fig.  (12)  presents  values  of  the 
computed  stresses  obtained  from  the  respective  analyses.  It  can  be  seen  that  near  the  region  where  the 
load  is  applied  (z  <  2.5  m.),  all  the  solutions  yield  the  similar  values  which  over  estimate  the  vertical 
stresses  indicating  effects  of  stress  concentration  are  dominant  within  this  region.  It  also  indicates  that 
the  coupled  boundary  elements  have  no  significant  effect  within  this  region.  The  relative  error  of  the 
stresses  are  shown  in  fig.  (13).  From  this  graph  we  can  see  that  the  prediction  of  the  stresses  in  the  far 
field  by  the  initial  finite  element  analysis  FEM  (100)  is  in  error  by  as  much  as  40%.  This  however  is 
reduced  to  less  that  4%  when  the  mesh  size  is  doubled  FEM  (169A).  The  coupled  solution  particularly 
the  unsymmetric  version  BEFEC  (100/18)  results  in  a  slightly  better  prediction  as  compared  with 
FEM  (169A)  where  stresses  are  within  2%  of  the  analytical  solution  at  the  far  field.  It  is  only  in  this 
particular  instance  that  there  is  a  observable  diflerence  between  the  symmetric  and  unsymmetric 
coupling  methods.  The  difference  is  significant  enough  such  that  the  symmetric  results  are  slightly  less 
accurate  than  the  those  obtained  using  the  expanded  mesh. 

Problem  2:  Uniformly  Distributed  Load  Applied  to  an  Elastic  Half-Plane.  To  solve  this  problem,  five 
different  analyses  were  considered  namely:  (1)  a  finite  element  analysis  consisting  of  a  40  m.  x  40  m. 
region  using  121  nodes  and  100  elements,  FEM  (121)  [see  fig.  (14)],  (2)  a  finite  element  analysis 
consisting  of  a  80  m.  x  80  m.  region  using  225  nodes  and  196  finite  elements,  FEM  (225A)  [see  fig. 
(15)],  (3)  a  finite  element  analysis  consisting  of  a  40  m.  x  40  m.  region  using  225  nodes  and  196  finite 
elements  FEM  (225B),  [see  fig.  (16)]  and  (4)  a  coupled  boundary  element/finite  element  analysis 
consisting  of  a  40  m.  x40  m.  region  with  121  nodes  and  100  elements  together  with  20  boundary 
elements  (BEFEC-121/20).  (5)  a  coupled  boundary  element/finite  element  analysis  similar  to  that 

performed  in  (4)  but  with  the  skew-symmetric  part  of  the  boundary  element  stiffness  matrix  discarded 
(BEFEC2-121/20).  Symmetry  in  geometry  was  considered.  The  material  properties  were  £=30  MPa, 
1/  =  0.25.  The  applied  load  had  a  magnitude  of  1  KN/m  and  was  applied  over  a  length  of  10  m. 

Results  were  compared  between  the  finite  element  solutions  where  the  far  field  displacements  were 
assumed  to  vanish  and  that  using  the  boundary  elements. 

The  following  observations  are  made. 

1.  Vertical  Displacements  along  the  centerline:  Because  of  symmetry  along  this  line  ,  the 
horizontal  displacements  are  zero.  Thus  only  vertical  displacements  were  considered.  The  results  arc 
presented  in  fig  (17).  They  show  that  using  the  conventional  method  of  analysis,  FEM  (121),  a 
significant  disparity  between  the  analytical  and  numerical  solution.  This  discrepancy  is  affected  very 


28 


(ed>|)  sssjis 


LO 

lO 

lO 

CO 

If) 

CVJ 

If) 

1— 

If) 

d 

d 

d 

CO 

d 

d 

CVJ 

d 

d 

d 

d 

o 

d 

29 


Fig.  12  Problem  1:  SIgma-YY  Along  Diagonal  Elements 
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Fig.  13  Problem  1:  Error  in  Sigma-YY  Along  Diagonal  Elements 
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(225B)  Mesh 


Analytical 


Y  Coordinate  (m.) 

Fig.  17  Problem  2:  Vertical  Displacements  Along  Centerline 


little  by  refining  the  mesh  as  was  done  in  FEM  (225B).  The  prediction  of  the  vertical  displacements 
along  the  centerline  improved  by  coupling  the  solution  with  boundary  element  method  BEFEC  ( 
121/20),  BEFEC2  (121/20)  or  increasing  the  dommn  of  the  problem  FEM  (225A)  the  latter  option 
giving  a  slightly  better  prediction  the  latter  giving  an  answer  very  close  to  the  analytical  solution.  As 
in  the  previous  problem  both  symmetric  and  unsymmetric  solutions  give  the  same  answer. 


2.  Horizontal  displacements  along  the  surface.  Evaluating  the  analytical  solution,  it  can  be  seen 
that  the  horizontal  displacements  along  the  traction  free  surface  are  a  constant  equal  to 
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From  fig.  (18),  it  can  be  noted  that  in  all  five  analyses,  there  is  a  substantial  discrepancy  between  the 
numerical  and  analytical  solution.  Because  in  the  finite  element  solutions,  the  horizontal  displacements 
are  assumed  to  vanish  away  from  the  loaded  area,  all  three  finite  element  solutions  diverge  from  the 
analytical  solution  as  pioints  away  from  the  loaded  area  are  considered.  Both  of  the  coupled  solutions 
give  displacements  which  are  close  to  the  analytical  solution  due  to  the  fact  that  the  above  property  of 
the  problem  is  modeled  using  the  fundamental  solution  for  the  boundary  element  method.  It  should  be 
stated  however  that  increasing  the  domain  of  the  problem  results  in  an  improvement  of  the  predicted 
horizontal  displacements. 


3.  Vertical  displacements  along  the  surface.  Results  presented  in  fig.  (19)  show  that  like  the 
displacements  along  the  centerline,  the  conventional  method  of  finite  element  analysis  FEM  (121) 
results  is  a  substantial  disparity  between  the  analytical  and  numerical  solution.  This  disparity  is  little 
affected  by  refining  the  mesh  as  was  done  in  FEM  (225B).  The  prediction  of  vertical  displacements  can 
however  be  improved  by  coupling  the  solution  with  the  boundary  element  method  BEFEC  (121/20),  or 
increasing  the  domain  of  the  problem  as  in  FEM  (225A).  In  this  particular  case,  the  displacements 
obtained  using  the  coupled  solution  are  much  more  accurate  than  the  finite  element  solution. 

4.  Horizontal  Normal  Stresses  (<t^^)  along  the  diagonal  elements.  From  fig  (20)  it  can  be  seen 
that  the  stresses  predicted  by  the  initial  finite  element  method  alone  FEM  (121)  is  slightly  over 
predicted  in  the  near  field  and  considerably  underpredicted  in  the  far  field.  It  is  noteworthy  to  point 
out  that  increasing  the  domain  as  in  FEM  (225A)  or  coupling  with  the  boundary  element  method 
BEFEC  (121/20)  increases  that  error  in  the  overprediction  of  the  stresses  near  the  loaded  area.  These 
error  can  be  attributed  to  the  approximate  representation  of  the  distributed  load  with  equivalent  nodal 
loads.  It  should  be  pointed  out  that  the  greatest  errors  in  this  particular  region  results  from  the 
coupling  of  the  finite  element  method  with  the  boundary  element  method.  In  comparison,  the  coupled 
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Fig.  18  Problem  2:  Horizontal  Displacements  Along  Surface 
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Fig.  19  Problem  2;  Vertical  Displacements  Along  Surface 
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Problem  2:  Sigma-XX  Along  Diagonal  Elements 


method  yielded  stresses  which  agreed  very  will  with  the  analytical  solution  in  the  far  field.  This  is  not 
true  for  the  values  obtained  using  the  finite  element  method  alone.  Significant  error  were  noted  in  the 
far  field  stresses.  From  the  relative  error  of  these  stresses  as  shown  in  fig.  (21)  it  can  be  seen  that  both 
analyses  result  in  error  of  up  to  40%.  It  should  also  be  emphasized  that  expanding  the  finite  element 
mesh  as  in  FEM  (225A)  does  not  improve  the  prediction  of  stresses. 

5.  Shear  Stresses  (o’^y)  along  the  diagonal  elements.  Results  of  the  analyses  performed  are 
presented  in  fig.  (22).  From  the  graph,  it  can  be  seen  that  there  is  a  significant  over  prediction  of  the 
stresses  in  the  region  close  to  the  loaded  area.  Also,  it  should  be  noted  that  all  method  resulted  in 
similar  predicted  stresses.  However,  in  the  far  field,  the  initial  finite  element  solution  FEM  (121) 
significantly  underpredicts  the  stresses.  This  underprediction  is  reduced  by  expanding  the  mesh  as 
shown  by  FEM  (255A).  Results  very  close  to  the  analytical  solution  are  obtained  when  the  solution  is 
coupled  with  the  boundary  element  method  as  in  BEFEC  (121/20)  and  BEFEC2  (121/20).  Comparing 
the  relative  error  as  presented  in  fig.  (23),  it  can  be  seen  that  the  FEM  (121)  far  field  solution  is  grossly 
in  error.  In  particular,  an  underprediction  of  60%  occurs  in  the  initial  finite  element  solution  FEM 
(121).  This  error  is  substantially  reduced  to  20%  by  expanding  the  mesh  to  twice  the  initial  size  as  in 
FEM  (225A).  The  results  obtained  using  both  the  unsymmetric  BEFEC  (121/20)  and  symmetric 
BEFEC2  (121/20)  coupled  solution  yield  identical  values  within  2%  of  the  analytical  solution. 

6.  Vertical  Normal  Stresses  (<ryy)  along  the  diagonal  elements.  Fig.  (24)  presents  values  of  the 
computed  stresses  obtained  from  the  respective  analyses.  It  can  be  seen  from  this  graph  that  the  effect 
of  the  approximation  of  the  distributed  load  with  equivalent  nodal  loads  results  in  a  significant 
underprediction  of  the  stresses  in  all  methods  used.  It  can  also  be  seen  from  this  graph  all  methods 
examined  yield  similar  values.  The  relative  error  of  the  stresses  are  shown  in  fig.  (25).  From  this  graph 
we  can  see  that  the  prediction  of  the  stresses  in  the  far  field  by  the  initial  finite  element  analysis  FEM 
(121)  is  in  error  by  as  much  as  30%.  This  however  is  reduced  to  less  that  4%  when  the  mesh  size  is 
doubled  FEM  (225A).  The  unsymmetric  coupled  solution  BEFEC  (121/20)  results  in  a  slightly  better 
prediction  as  compared  with  FEM  (225A)  where  stresses  are  within  3.5%  of  the  analytical  solution  at 
the  far  field.  It  can  also  be  seen  that  there  is  a  significant  difference  in  the  symmetric  and  unsymmetric 
coupling  technique  with  respect  to  the  stresses  computed  in  the  far  field.  In  the  symmetric  solution,  the 
stresses  in  the  far  field  were  in  error  by  as  much  as  5%.  It  can  be  concluded  that  in  the  prediction  of 
the  vertical  normal  stresses,  expanding  the  mesh  and  symmetric  coupling  of  the  boundary  element 
method  result  in  similar  results. 
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Fig.  21  Problem  2:  Error  in  SIgma-XX  Along  Diagonal  Elements 
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Fig.  22  Problem  2:  Sigma-XY  Along  Diagonal  Elements 
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Fig.  23  Problem  2:  Error  in  Sigma-XY  Along  Diagonal  Elements 
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Fig.  24  Problem  2:  Sigma-YY  Along  Diagonal  Elements 
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Fig.  25  Problem  2:  Error  in  Sigma-YY  Along  Diagonal  Elements 


Conclusions  and  Recommendations 


The  following  observations  can  be  made  with  respect  to  the  results  of  the  proposed  method  of 
analysis. 

1.  The  coupling  of  the  boundary  and  finite  element  method  provides  sm  effective  method  for 
solving  soil-structure  interaction  problems.  The  method  has  been  shown  to  be  efficient  in  solving  linear 
problems  involving  infinite  half-plane  domains.  Both  the  displacements  and  stresses  obtained  using  this 
method  are  in  close  agreement  with  existing  analytical  solutions.  It  should  be  mentioned  that  the 
choice  of  the  fundamental  solution  used  in  the  boundary  element  method  accounts  for  the  efficiency  of 
the  method.  This  is  because  ail  degrees  of  freedom  in  the  boundary  element  system  are  assembled  in 
the  finite  element  system.  Thus  the  size  of  the  boundary  element  system  is  kept  to  a  minimum. 

2.  It  was  noted  that  for  the  problems  solved,  the  difference  between  using  the  symmetric  and 
unsymmetrir  part  of  the  boundary  element  stiffness  matrix  is  almost  negligible  with  respect  to  the 
prediction  of  displacements.  Thus  it  is  possible  to  discard  the  skew-symmetric  portion  of  the  boundary 
element  stiffness  matrix  and  still  obtain  accurate  answers. 

3.  Increasing  the  domain  of  the  mesh  does  significantly  reduce  the  error  in  the  prediction  of  the 
far  field  vertical  normal  and  shear  stress  but  does  not  reduce  the  error  in  the  prediction  of  the  far  field 
horizontal  shear  stress.  In  comparison,  the  use  of  the  coupled  method  results  in  reduced  error  for  all 
components  of  stress.  Additionaly,  these  errors  are  significantly  less  than  those  resulting  from  the  finite 
element  method  alone. 

4.  It  should  be  emphasized  that  in  the  cases  where  identical  results  are  obtained  by  expanding  the 
finite  element  mesh  and  coupling  the  finite  element  method  with  the  boundary  element  method,  the 
latter  approach  is  still  more  efficient  because  it  results  in  far  more  less  equations  to  be  solved.  It  should 
also  be  pointed  out  that  the  results  obtained  by  symmetric  coupling  are  at  worst  as  accurate  as  those 
obtained  by  expanding  the  mesh  and  at  most  times  are  as  accurate  as  the  unsymmetric  coupling 
technique. 

5.  In  the  problem  involving  the  concentrated  load,  it  was  observed  that  the  stresses  obtained 
using  the  finite  element  method  was  sufficiently  accurate  and  the  only  significant  errors  occurred  in  the 
far  field. 

6.  There  is  a  significant  discrepancy  which  results  from  the  representation  of  the  distributed  load 
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with  equivalent  nodal  loads.  This  discrepancy  cannot  be  removed  by  either  coupling  the  solution  with 
the  boundary  element  method  or  by  modifying  the  mesh  and  these  errors  can  even  be  intensified  by 
using  such  methods. 
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APPENDIX  A 

Analytical  Solutions  to  E^xample  Problems 


In  this  appendix,  the  complete  analytical  solution  to  the  problems  solved  are  presented.  These 
consist  of:  (1)  A  vertical  concentrated  load  applied  to  the  surface  of  jui  elastic  half-space,  (2)  A 
distributed  load  applied  to  the  surface  of  an  elastic  half-plane. 

A.l  Concentrated  Load  Applied  to  an  Elastic  Half-Plane 

The  displacements  and  stresses  corresponding  to  this  problem  are  given  by 


2  H  ^xx "  2  ^xx  ~~  2 

where  the  load  Py  is  applied  at  the  origin,  while  the  displacements  and  stresses  are  evaluated  at  the 
point  (x,y)  and  r  is  the  distance  from  the  origin  to  the  point  (x,j/)  such  that  r  =  .^i2  +  y2.  The 
constant  G  is  the  shear  modulus,  u  is  poisson’s  ratio  and  L  is  an  arbitrary  constant  chosen  so  that  at 
the  point  (L,0),  =  0. 

A. 2  Distributed  Load  Applied  to  the  Surface  of  an  Elastic  Half-Plane 

The  displacements  and  stresses  corresponding  to  this  problem  are  given  by 

tij  =  [(1  -  2j/){(x  -  a)  -  (x-f-o)  e2-ira)}  +(l-i/)y  ln(rj/r^)] 

“v  =  ^  [-(1  -2j/)  y  («j  -^2)  +(!-*')  {(*-“)  !«(»•?)-(*  +  «) 

-f-  (L  a)  In  (L  -t-  o)^  —  (L  —  o)  In  {L  —  o)^}] 
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^  [^1  -  ^2  +  y(*  -  ")/'■!  ■" 

^VW  =  -^  [^1  -  ^2  -  V(*  -  O)/'-?  +  V(*  +  “)/’‘2] 

-x„  =  ^[lA?-Vr^2] 

In  these  expressions,  the  following  abbreviations  are  used 
=  arctan  a  ^2  =  arctan 

rj  =  (x-a)2  +  y^  =  (t+a)^  +  y^ 

where  the  distributed  load  w  is  applied  from  points  (— a,0)  to  {a,0),  while  the  displacements  and 
stresses  are  evaluated  at  the  point  {x,y)  and  r  is  the  distance  from  the  origin  to  the  point  (i,y)  such 
that  r  =  ^x2  +  y2.  The  constant  G  is  the  shear  modulus,  v  is  poisson’s  ratio  and  X  is  an  arbitrary 
constant  chosen  so  that  at  the  point  (X.O),  =  0. 
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APPENDIX  B 

Fundamental  Solution 


In  this  app>endix  the  fundamental  solution  for  the  displacements  and  stresses  due  to  a  unit  load 
applied  to  a  half-plane  is  presented.  The  complete  solution  consists  of  two  parts  namely:  (1)  the  Kelvin 
solution  representing  the  solution  to  an  infinite  plane,  (2)  the  complementary  solution  which  when 
added  to  the  Kelvin  solution  gives  the  half-plane  solution.  Thus  the  complete  solution  can  be  written  as 
follows: 


_ L_ _ 

8tG(1-i/)  [r2 


““  "  8»C('l-,-)  l"(r)  +  ^ 
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Figure  B.1  Unit  Load  Appiied  to  an  Eiastic  Haif-Piane 


-  1 
4t(1  —  i/)r 


^^{(1-2.)®  +  ^} 


fi^{(l-2.)(f)  +  ^} 


J^{(1-2.)(S)  +  ^} 


B.2  Complementary  Solution 


[8(1-i/)2  -{Z-4u)]ln  R  + 


[(3  — 4»/)iiJj  —2xx]  4x  X 

- F - *~R^ 


-  ,,  f  (3  — 4i/)r,r,  4z  x  1 

4(1 ->/)(! -2./)g I 


r  f  (3  — 4»/)r,r5  4x  x  iZ.r,  ^  ^  1 

«!,  =  - +  ■1(1  —Ki  -2-)« 


J  - (8(1  - .f  - (3 - 4.))  I„  «  +  I(3-4.H-2.?1^4^ 


’5„  =  -  K,  { 


{3x  +  x)(l  -  2|/)  ,  +  2x  x)-(4xr2)(l -2|/)  16x  x  iJjTj 

_o  ^  “  _c 
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/  (l-2i/)  2[  *2  _2ix  -x2+  2xRi(l-2i/)]  16r  x  ie?\ 

_  - - +  — F“/ 


=  -K  I 

I  R 


3x)(l-2z/)  2[  R^(rl  +  2x^)  -  2xr^  +  2xr^(l  -  2t/)]  16x  x  R] 


—  _  tr  ^  j  ~  2*^)  2[  x^  —  x^  +6xx  —  2xftj(l  —  2i/)]  16x  x  r2\ 

'  R* 

^  -  IT  /  (3?  +  a:)(l  -  2i/)  ,  2[(2x  x  +  -  2xR^(l  -  2v)]  16x  i 

'1  B?  +  ^6 - 1 


where 


^%iS  -  ~  Bjr.^  I 


3(1  -2t/)  _  2[  rl  -  4xx  -  2x^-  2xi?,(l  -  2t/)]  16i  x  Rf 


0  =  arctan 


(«.) 


_  _ 1 _ 

^  [8x(l-j/)G] 


_  1 

•  '  [4x( !-«.)] 
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APPENDIX  C 

Derivation  of  Transformation  Matrix 


From  the  finite  element  method,  the  equivalent  nodal  loads  f-j  for  a  given  traction  distribution 
h((,  Tf)  over  the  boundary  F  are  given  by 


/.,=  /  dr  (C.l) 

J  p 

where  is  the  ith  traction  compoenent,  is  the  lagrangian  interpolating  functions 

corresponding  to  the  jth  node.  Because  of  the  properties  of  the  interpolating  function,  only  to  those 
components  corresponding  to  nodes  on  the  boundary  F  are  non-zero.  Furthermore,  these  non-vanishing 
comp>onents  degenerate  to  functions  in  one  variable.  Thus,  the  above  integral  can  be  rewritten  as 


/.,=  /  hiiO^jiO  dr 
•'  r 

where  ^  is  the  coordinate  along  the  boundary  F. 


(C.2) 


Similarly,  from  the  boundary  element  method  it  is  assumed  that  the  tractions  over  the  entire 
element  vary  in  the  following  manner 

MC)=  V,(C)  (C-3) 

where  h  -  is  the  ith  component  of  the  traction  at  the  jth  node. 

Thus,  combining  the  two  above  equations,  the  equivalent  nodal  loads  can  be  obtained  by  the 
cition 

/,^  =  h..,7V^^  (C.4) 

where 

/  MO^jiOdC  (C.5) 

J  p 

For  the  coupling  between  a  bilinear  quadrilateral  finite  element  and  linear  boundary  element,  the 
corresponding  shape  functions  are  given  by 


V'i(C)  =  |(l-C)  t^2(C)  =  5(H-0 


(C.6) 
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If  the  boundary  T  is  a  straight  line  of  length  /,  the  integral  in  the  above  equal  evaluates  to 


(C.7) 


For  a  general  boundary  element  system,  the  above  matrix  must  be  assembled  into  the  global  N 
matrix  very  similar  to  the  way  line  elements  are  assembled  in  a  finite  element  system. 
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